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Abstract 


The computationally efficient HZETRN code has been used in recent trade studies for lunar and 
Martian exploration and is currently being used in the engineering development of the next 
generation of space vehicles, habitats, and extra vehicular activity equipment. A new version 
(3DHZETRN) capable of transporting High charge (Z) and Energy (HZE) and light ions 
(including neutrons) under space-like boundary conditions with enhanced neutron and light ion 
propagation is under development. In the present report, new algorithms for light ion and neutron 
propagation with well-defined convergence criteria in 3D objects is developed and tested against 
Monte Carlo simulations to verify the solution methodology. The code will be available through 
the software system, OLTAR1S, for shield design and validation and provides a basis for personal 
computer software capable of space shield analysis and optimization. 

Introduction 

High charge ( Z) and Energy (HZE) ion transport developed parallel to the experimental studies of 
Schimmerling et al. [1986] at the Lawrence Berkeley Laboratory and involved solving the Boltzmann transport 
equation for mono-energetic ion beams in the context of the continuous slowing down approximation (CSDA) 
[Wilson et al. 1984]. The main purpose was to identify the source of computational limitations. This was found in 
part to be the inadequacy of available nuclear data through comparing computational results with the experimental 
ionization data for a broad beam of 20 Ne ions [Wilson et al. 1984]. Emphasis was soon overtaken by the need to 
establish the scope of the galactic cosmic ray (GCR) protection problem, and marching procedures were used to get 
first order estimates of shielding requirements [Wilson and Lamkin 1975, Lamkin 1974, Wilson and Badavi 1986, 
Wilson et al. 1991]. Testing the new computational marching model against atmospheric air shower data again 
pointed to the inadequacy of the available nuclear data [Wilson and Badavi 1986], and development of a semi- 
empirical nuclear fragmentation model (NUCFRG) followed, leading to improved comparisons with experimental 
deep penetration data [Wilson et al. 1987a, 1987b]. This new NUCFRG model likewise improved comparisons with 
Schimmerling’ s laboratory studies [Shavers et al. 1993, Wilson et al. 1991]. The next dozen years emphasized 
spaceflight validation of the marching solution and nuclear fragmentation model improvements [Wilson et al. 1995, 
1998; Badhwar et al. 1996; Shinn et al. 1998; Cucinotta et al. 1998; Hugger et al. 2003; Nealy et al. 2006; Wilson et 
al. 2005, 2006, 2007; Slaba et al. 201 la, 2013; Norman et al. 2012, 2013]. 

Advanced solution methods of the Boltzmann equation continued to develop [Wilson et al. 1994a, 1994b] 
but only slowly after NASA interest in deterministic transport code and related nuclear model development shifted 
in favor of Monte Carlo methods [Armstrong and Colburn 2001, Pinsky et al. 2001]. Focus of the Langley Research 
Center activity turned to trade studies, flight validation, and development of engineering design methods including 
shield optimization processes for which deterministic methods are very successful [Wilson et al. 2004a]. Renewed 
interest within NASA for deterministic code development is giving new emphasis for improved solution methods 
[Wilson et al. 2003a, 2006; Tweed et al. 2006a, 2006b; Slaba et al. 2010a, 2010b; Norman et al. 2013] with a focus 
on nuclear modeling activity as well [Cucinotta et al. 1994, 1998, 2007; Norbury et al. 2007; Adamczyk et al. 2012]. 
As a result, recent developments mainly utilized the older semi-empirical NUCFRG2 model with revisions [Wilson 
et al. 1994c, 1995, 2006] including energy downshifts and momentum dispersion [Tripathi et al. 1994]. Comparison 
of NUCFRG2 with other models and nuclear fragmentation experiments are given by Zeitlin et al. [1997, 2011], 
Golovchenko et al. [1999], Cucinotta et al. [2007], and Walker et al. [2006]. The present development will use the 
more recent NUCFRG3 [Adamcyzk et al. 2012] with an added light ion (Z < 2) coalescence model and additional 
electromagnetic dissociation channels. 

To support the need for keeping records of astronaut exposures in space operations, considerable time has 
been spent in validating and testing in Shuttle flight operations with an array of instruments [Badhwar et al. 1995, 
1996, 2001]. As a result of the need to design and test new systems for future exploration, there has been renewed 
interest in flight validation using the International Space Station (ISS) as a measurements platform in low Earth orbit 
(LEO) [Hugger et al. 2003; Wilson et al. 2005, 2006, 2007; Nealy et al. 2006; Slaba et al. 201 la, 2013]. 

In the present report, the current status of transport code development is reviewed with emphasis on 
extending these developments into a 3D version of the HZETRN code. These advances will use available Monte 
Carlo codes Geant4 [Agostinelli et al. 2003, Geant4 Collaboration 2012a, 2012b], FLUKA [Fasso et al. 2005, 
Battistoni et al. 2007], and PHITS [Sato et al. 2013] to judge the veracity of these developments especially with 
regard to their 3D aspects. Upon completion, the 3DHZETRN code will be integrated into the web based OLTARIS 
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[Singleterry et al. 2011] software system for general testing, validation against spaceflight data, and use in spacecraft 
design. 

Deterministic Code Development 

The relevant transport equations are the time-independent linear Boltzmann equations derived on the basis 
of conservation principles [Wilson et al. 1991] for the flux density ^ (x, IT, E ) of type j particles as 


cr jk (E,E',£l,n, ')4> k (x,Sl\E')dfi'dE' - ^ (£)</>. (a, ft, £), 

k 

and solved subject to a boundary condition over the enclosure of the solution domain as shown in Fig. 1. 


( 1 ) 



Fig. 1. Geometric relations of quantities useful in solving equation (1). The symbol n is a unit normal vector. 

The Oj(E) and o jk (E,E\£l,£l') are the media macroscopic cross sections. The double differential cross section, 
o jk (E,E ’,£ 2 , 12 '), represents all those processes by which type k particles moving in direction £ 2 ' with energy E ' 
produce a type j particle in direction £2 with energy E (including decay processes). Note that there may be several 
reactions that produce a particular product, and the appropriate cross sections for equation (1) are the inclusive ones. 
Exclusive processes are functions of the particle fields and may be included once the particle fields are known. The 
total cross section oj(E) for the medium may be written for each particle type as 

a j (E) = <xj a<) (E) + af l) (E) + af uc) (E) , (2) 

where the first term (at) refers to collision with atomic electrons, the second term (el) is for Coulomb elastic 
scattering on the atomic nuclei, and the third term (nuc) describes nuclear scattering and reactions where the minor 
nuclear inelastic processes (excitation) have been ignored. The corresponding differential cross sections are 
similarly ordered. Many atomic collisions (~10 6 ) occur in a centimeter of ordinary matter, whereas ~10 3 nuclear 
Coulomb elastic collisions occur per centimeter, while nuclear processes are separated by a fraction to many 
centimeters depending on energy and particle type. Solution methods for equation (1) use perturbations based on 
the ordering of the physical cross sections with atomic interactions as the first physical perturbation; special methods 
are used for neutrons for which atomic cross sections are zero. Coulomb scattering is the second order physical 
perturbation process. The nuclear processes are treated in detail as the third order physical perturbation. 

The atomic/molecular interaction terms of equations (1) and (2) may be written as [Wilson et al. 1991] 
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( 3 ) 


E f 4 yf(E,E',n,Cl')ct> k (x,ri',E')dn'dE' - ^(Etyfa^E) 

= E^n( E + A-h^ix^E + A~h n ) - af\E)^(xM,E) 


E 

n l 

4-i d 




< 7 <* l (£)^(*,n,i ?y 


-^"»(W j (x,n,£) + o(^) 




where ^ is the mass of particle type y, and the stopping power S/E) is given by 

S j (E) = E s n^(E)- 

n 


(4) 


In equations (3) and (4), n denotes the electronic excitation levels, and s n refers to the corresponding excitation 
energies. Evaluation of the stopping power by equation (4) appears deceptively simple. However, all of the excited 
states including continuum states of the atomic/molecular system need to be evaluated. Furthermore, the projectile 
remains a bare ion except at low energies where the projectile ion atomic orbital states begin to resonate with the 
relative speed of the electrons of the media leading to electron capture and lowering of the ion charge (the stopping 
power is then related to RMS charge, called effective charge). These issues are further discussed in Wilson et al. 
[1991, 2005, 2006] and Tai et al. [1997]. The higher order terms of equation (3) are neglected in the CSDA and are 
further discussed elsewhere in connection to mono-energetic ion beams [Wilson et al. 2005, 2006]. The second 
physical perturbation, Coulomb scattering, results in angular diffusion about the ion direction leading to small 
angular corrections that are of little interest in the nearly omnidirectional space environment [Wilson et al. 1974, 
2005, 2006] and are therefore ignored. Henceforth, the total cross sections and the double differential cross sections 
appearing in equation (1) are the nuclear scattering and reaction values, and the superscript ( nuc ) will not be written. 

Equation (1) can be written in the CSDA, ignoring multiple Coulomb scattering, as 


B [<t>.(x,n,E)} = zi;l a jk (E,E',n,n')^ k (x,n',E')dn'dE' , (5) 

k 

where the differential operator on the left hand side is defined as 

b [0 . ( x , n ,E)] = n. v^. (*, n, e) - 2- A [ s . ( e)^ (*, «,£)] + ( x , n, e) . (6) 

The right-hand side of equation (5) excludes the atomic/molecular processes now appearing on the left as the energy 
shifting operator in addition to the spatial drift and nuclear absorption terms. Neutral particles would have null 
atomic cross sections (except for negligible magnetic coupling) for which the stopping power term of B in equation 
(6) is set to zero. Application of CSDA in both laboratory and space shielding has been wide and the resulting errors 
discussed elsewhere [Tweed 2006a, 2006b]. Equation (5) can be rewritten as a Fredholm equation [Wilson 1977] 
given by 


Ax,n,E) = 




Sj(E)Pj(E) 


[r (n,*),n,E 


E c /“ f 4n p 3 ( E >jki.E\ E ", o, n ')4> k (x n ,n\ E ")dn 


(7) 


SAE)PAE)^Je J E' Jin 


y dE"dE\ 


with 
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E i = R~ 1 \P - d + R^E)}, 


( 8 ) 


Xsi =x + [R.(E) - RjiE'W . 


( 9 ) 


In equations (7) and (8), T(Q,a:) is the point on the boundary connected to x along -Q, p is the projection of jc onto 
12 (Fig. 1), d is the projection of r(Q,x) onto Q, Rj(E) is the distance a type j ion with energy E will travel before 
losing all of its energy to excitation/ionization of atomic electrons, and Pj(E) is the probability that a type j ion with 
energy E will not have a nuclear interaction in coming to rest in the media. The zero charge limit of equation (7) is 
discussed elsewhere [Wilson et al. 1991]. One obstacle to solving either equations (5) or (7) is the need to evaluate 
the integral d£l' at arbitrary locations within the media and development of methods to handle this limitation 
expediently in ways amenable to efficient computational practice. 

The range-energy relationship is given (Fig. 2) by 


R ,(E) = Aj 


E dE ' 


S+E' 


( 10 ) 


and the nuclear attenuation function is given (Fig. 3) by 


Pj{E) = exp 


a p -JEl 

J d 0 Sj(E') 


dE' , 


( 11 ) 


where oj(E ? ) includes nuclear elastic and reactive processes. 

The nuclear reactive and scattering differential cross sections can be written in the following form 


a jk (E,E',n,n') = ±-a jk .J E ,E') + a jkfor (E,E',n,n'), ( 12 ) 

where the first term (iso) is isotropic and associated with lower energy particles produced, including target 
fragments, and the second term (for) is highly peaked in the forward direction (Q ~ Q f ) and is associated mainly 
with direct quasi-elastic events and projectile fragmentation products [Wilson 1977, Wilson et al. 1988]. 



Fig. 2. Range-energy relationship of ions in aluminum. Results for Z = 1 and Z = 2 are nearly overlapping. 
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Fig. 3. Probability of nuclear reaction of ions in aluminum. Results for Z= 1 and Z = 6 are nearly overlapping. 



Fig. 4. Isotropic (iso) and forward (for) neutron spectra produced by 500 MeV protons in aluminum. 


Surprisingly, even nucleon-induced reactions follow this simple form, and the isotropic term extends to relatively 
high energies (see Fig. 4). For nucleon-induced reactions, the following form of the differential cross section has 
been used in versions of FLUKA [Ranft 1980] as follows 

cr jk (E,E'MM') = v^E'^E'V^E'^E,^), (13) 


where Vjk(E') is average multiplicity, fjk(E,E ’) is the energy spectrum of produced fragments, and the Ranft angular 
factor used in early versions of FLUKA is 


E,Ap) 


N r exp(— 0 2 / X R ) ,O<0<7r/2 
N r exp(- 7 r 2 / 4A^) , tt / 2 < 0 < tt 9 


(14) 


with cos 6 = £l £l f and 
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X R = (0.12 + 0.00036^) / E . 


(15) 


In equations (14) and (15), A T is the mass of the target nucleus, and Nr is a normalization constant. The Ranft 
formalism was modified by adding a quasi-elastic component [Wilson et al. 1988] and was instrumental in achieving 
the correct neutron KERMA [Wilson et al. 1991]. The isotropic component of equation (12) is determined by 
integrating equation (13) over the back hemisphere (i.e. cos(0) < 0) according to 

"friso&E') = 2 f 2 ^ B a jk (E,E',n,n')dn', (16) 


where 2nB refers to the back hemisphere. The separation in phase space given by equations (12) - (16) is exploited 
in obtaining corrections to the straight-ahead approximation for light particles (Z < 2). The heavy ion (Z > 2) 
projectile fragment cross sections are further represented by 


* jk (E,E',n,n') = N t 


^f2exp[-2mV^(l - fi.fi') / s t . k ]exp[-(£ - E'+ X jk ) 2 / 2 e% 

P n£ jk 


(17) 


where m is the nucleon rest mass, A Jk is related to the momentum downshift, s jk is related to the longitudinal 
momentum width, s t>jk is related to the transverse momentum width, N t is the transverse normalizing factor, and the 
fragmentation cross sections, o jk (E ' ) are obtained herein from NUCFRG3. Since the transverse width is small 
compared to the projectile-fragment energy, the transverse function is highly peaked about the forward direction 
(Q-lT-l). 

Atomic interactions limit the contributions of charged particles in the transport process. For example, the 
protons and alpha particles produced in aluminum below 100 A MeV contribute to the fluence only within a few 
centimeters of their collisional source and the heavier ions are even more restricted (see Fig. 2). This is an important 
factor in that the transported secondary charged particle flux tends to be small at low energies and the role of 
additional nuclear reactions by these low-energy ions are likewise limited (see Fig. 3), an important fact to be 
implemented later in the current development. 

Straight-ahead solution (N = 1) 

The approach to a practical solution of equation (5) is to develop a progression of solutions from the simple 
to increasingly complex, allowing early implementation of high-performance computational procedures and 
establishing a converging sequence of approximations with established accuracy criteria and means of 
verification/validation. The first step leading to the lowest order solution uses the straight-ahead approximation as 
guided by the nucleon transport studies of Alsmiller et al. [1965] using Monte Carlo methods in which the 
differential cross sections were approximated as 

<T. k (E,E\n,Sl') = <T. k {E,E')6{fl-Sl') 9 (18) 


resulting in dose and dose equivalent per unit fluence to be within the statistical uncertainty of the Monte Carlo 
result obtained using the fully angle dependent cross sections in slab geometry. The relation of angular dependent 
cross sections to spacecraft geometry in space applications was further examined by Wilson and Kandelwal [1974] 
using asymptotic expansions about angular divergence parameters, demonstrating errors in the straight-ahead 
approximation to be on the order of the square of the ratio of distance of divergence (few centimeters) to radius of 
curvature of the shield (few to several m), resulting in a small error ~10" 4 in most human rated space systems 
[Wilson and Khandelwal 1974]. Introducing approximation (18) into equation (7) reduces the Fredholm equation to 
a Volterra equation that can be solved using marching procedures [Wilson and Famkin 1975, Wilson 1977, Wilson 
and Badavi 1986]. Considering the success of the straight-ahead approximation for nucleons encourages its use in 
HZE transport where the mass of heavier ions reduces even further the angle of scatter. The success of equation (18) 
results from the fact that lateral diffusion from a given ray is fully compensated in a flat plate (i.e. slab with infinite 
lateral dimensions) by lateral diffusion along adjacent rays. The verification and validation processes are described 
elsewhere [Wilson et al. 2005, 2006]. 
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Initial quantification of the usefulness of the code evaluated dose and dose equivalent (as given by the 
ICRP 26 [ICRP 1977]) in 30 cm slab of water behind a 20 g/cm 2 slab shield of aluminum (and alternately of iron) 
for the Webber [1966] approximation of the February 23, 1956 solar particle event (SPE) spectrum given as 


<t> P W 


IQ 9 )# + 938) 
200 p{E) 


exp{[239.1 — p(.E')] / 100} , 


(19) 


with p(E) = -\Je(E + 1876) . Comparisons were made against 3D Monte Carlo results (an early version of HETC) 

by Scott and Alsmiller [1967]. The results are shown in Fig. 5 for dose and Fig. 6 for dose equivalent. Statistical 
uncertainty of the Monte Carlo result is on the order of ten to fifteen percent, and the early version of HZETRN 
certainly lies within these bounds. 




Fig. 5. Dose in water target protected by 20 g/cm 2 aluminum shield (left pane) and iron shield (right pane) exposed 

to Webber SPE spectrum (from Wilson et al. [1991]). 



Fig. 6. Dose equivalent [ICRP 1977] in water target protected by 20 g/cm 2 aluminum shield (left pane) and iron 
shield (right pane) exposed to Webber SPE spectrum (from Wilson et al. [1991]). 
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The newly developed dynamic LEO environmental models, shield models, and transport algorithm are 
shown to have an end-to-end uncertainty estimated to be 20 percent by Badhwar et al. [1995, 1996, 2001]. Aside 
from its accuracy, this version of HZETRN provides a computationally efficient procedure allowing rapid 
dosimetric evaluations and providing a basis for development of advanced engineering design methods capable of 
design optimization processes [Wilson et al. 2003a, 2004b]. A sample of results is given in Table 1 and Fig. 7. A 
sequence of corrections to the straight-ahead approximation is now taken up for light ions and neutrons. The heavy 
ion nuclear cross sections, represented by equation (17), indicate that the straight-ahead approximation is suitable for 
heavy ions for most applications, including space, as evidenced by laboratory validation with mono-energetic beams 
[Shavers et al. 1993, Shinn et al. 1998, Walker et al. 2006]. 



Fig. 7. Observed dose rate versus calculated dose rate [Badhwar, Chapter 21 of Wilson et al. 1997]. 


Table 1. Compari son of straight-ahead HZETRN transport model with Shuttle flight data [Wilson et al. 2003b]. 
Flight Date DRNM* DLOC TLDf (pGy/d) Calc. (pGy/d) 


STS-41A 

11/83 

6421 

3 

64.6 

59.6 

STS-51D 

4/85 

6661 

4 

917.4 

889.3 

STS-31 

4/90 

5701 

1 

2141 

2290 

STS-43 

8/91 

5894 

4 

20.7 

18.6 

STS-62 

3/94 

6771 

1 

94.3 

89.2 

STS-65 

7/94 

6822 

2 

28.3 

25.1 

STS-67 

3/95 

6925 

3 

250.8 

238.1 

STS-80 

11/96 

6973 

4 

264.4 

256.5 

STS-82 

2/97 

7074 

1 

2978 

3080 

STS-91 

6/98 

6894 

1 

89.1 

83.2 

STS-101 

5/00 

6460 

2 

140.8 

131.1 

STS-92 

10/00 

6417 

2 

165.9 

153.4 


* Deep River Neutron Monitor count rate, | GCR corrected TLD-100 data 


Bi-directional neutron transport (N = 2) 

It is clear that the straight-ahead approximation provides accurate dosimetric results and a good 
approximation to the particle fluence in many circumstances (see Badhwar et al. [1995]). One limitation of the 
straight-ahead approximation is near the front boundary on a flat plate where the calculated neutron fluence vanishes 
(unless they are part of the external radiation environment). This results from the straight-ahead approximation itself 
and can be improved by connecting to 3D corrections containing backward leakage to the front boundary [Marchuk 
and Lebedev 1986]. This limitation is resolved by treating the leakage that is an important process of neutron 
transport in space systems but of less importance to the charged particle transport that is limited by range/energy 


8 


relations. This bi-directional approximation provides the first correction to the straight-ahead approximation and 
provides major improvements to the neutron fluence spectra. 

In the early transport studies of Wilson and Lamkin [1975] and Lamkin [1974], it was demonstrated that 
once neutrons are produced by a proton beam, the recoupling of neutrons back to the proton field is limited, but the 
neutron fields tend to build up with increasing penetration depth. This mainly occurs since many protons produced 
by neutron induced reactions are of lower energy than the forward components of equation (12) and are quickly 
removed by atomic interactions (especially these isotropically produced protons). The probability of even a 100 
MeV proton to undergo a nuclear reaction is only several percent (Fig. 3). The low energy neutrons on the other 
hand have no effective atomic interactions and propagate until a nuclear interaction occurs or they escape from the 
media through a boundary. Yet, these low energy protons and other light ions produced in tissue through neutron 
collisions contribute significantly to biological injury and must be included. This was the basis of studies starting 
with Clowdsley et al. [2000, 2001] and Heinbockel et al. [2003], with recent improvements provided by Slaba et al. 
[2010b]. Herein, we will concentrate on its extension to 3D transport methodologies. Future work will be a direct 
evaluation of equation (7) quantifying the propagated error of the 3D neutron code. 

A bi-directional neutron transport code was first investigated by Clowdsley et al. [2000, 2001] using 
multigroup methods and recoupling to the light ions using an analytic solution as a quadrature over the light ion 
collision source term. These types of approximations were further studied and improved by Slaba et al. [2010b] and 
later compared to Monte Carlo simulations in slab geometry [Heinbockel et al. 2011a, 2011b] and on the lunar 
surface [Slaba et al. 2011b]. The model has been implemented into OLTARIS and is the basis of the ray-by-ray 
transport methodology [Slaba et al. 2011b, Singleterry et al. 2011]. In that transport model, the charged particle 
components were first solved in the straight-ahead approximation with the neutron coupling represented only by the 
forward component of the differential production cross section. The isotropically produced neutrons are then solved 
using a coupled bi-direction (forward/backward) transport to complete the neutron solution. These isotropically 
produced neutrons further collide with the media, giving rise to an isotropic source of ions that are transported 
assuming only atomic interactions as an approximation to the complete solution. Neglected in this solution is the 
minor coupling between low energy charged particles and any produced fragments. This assumption is generally 
accurate since the range of the low energy isotropic ions is much smaller than their nuclear mean free path lengths 
(see Figs. 2 and 3). 

The governing equation in this approach is first given by equation (5). The solution methodology requires 
particle fluxes to be separated into forward (^ 7 / or ) and isotropic components (<^ 0 ) according to 

4>j (x, n, E) = 4> jJor (x, n, e) + 4> jMo (*, n, e) . (20) 

The differential cross section for neutron fragments (j = n) is also separated into forward and isotropic components 
according to equations (12) and (16). The differential cross section for charged fragments is treated in the straight- 
ahead approximation by equation (18). Upon substitution of the separated differential cross sections and particle 
fluxes into equation (5), the forward component of the charged particle (j ^ n) fluxes is obtained by solving the 
straight-ahead equation 


B[cj) jJor (x,n,E)] = f™a jk (E,E')<fi kJor (x,n,E')dE' , (21) 

k 

and the forward neutron component is obtained by solving 

B [<t> nJbr (x,Sl,E)] = J2 f™* nkJor (E,E')<t> kJor (xME')dE', (22) 

k 

for which the stopping power term in the operator B is zero for neutrons. Note that the full cross section of equation 
(18) appears in the ion transport of equation (21), but only the forward component of the cross sections appears in 
the transport of neutrons in equation (22). Equations (21) and (22) are solved simultaneously using methods 
previously developed [Wilson et al. 2006, Slaba et al. 2010a]. 

The forward fluxes give rise to an isotropic source of low energy neutrons which are treated by solving the 

equation 
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( 23 ) 


B Kieo&n> E )\ = f^f 4 ^nn^ E '^')^Jx^',E')dn'dE' 

+£ f E Jl °nk,iso (E, E',n,n ') 4 > kJor {x, Sl\E')d£l'dE'. 

k 

Whereas equations (21) and (22) take on value only if an external source is present at the boundary, equation (23) 
has no external sources, and the resulting field develops from the internal isotropic neutron source induced by the 
forward solution. Development of an algorithm for solving equation (23) is expedited by separating the isotropic 


neutron fluence into forward and backward propagating components according to 

Kso (*> E) = 4> n . so{f) (a, n, E) + <t> n . so{b) (x,-n,E), (24) 

and approximating the double differential cross sections as 

<7 nn (E,E',Q,Sl') = a nnJ (E,E')6(n - fi') + a nnJ) (E, E ')«(« + O'), (25) 

°nnA E ,E') = f 2vF a nn (E,E',n,n')dn' = a nnJor (E,E') + 2™ nn .jE,E') / 4,7, (26) 

°nn, b (E,E') = f 2vB <T nn (E,E',n,n')<m' = 2™ nn . S0 {E,E') / 47 r, (27) 

^KMf)(E,E') = f 2iTF a nk . S0 (E,E’,n,n’)dn' = 2na nk .jE,E') / 4tt, (28) 

VnkM^E’E') = f 2 ^n k ,,JE,E',n,n')dn' = 2ira nk .jE,E') / 4n . (29) 


In equations (26) - (29), 2nF and 27 l 5 represent integration over the forward and backward hemispheres with respect 
to H', respectively. The subscript if) and (b) in equation (24) refers to the forward hemisphere and backward 
hemisphere components of the isotropic neutron flux, respectively. Using equations (25) - (29) in equation (23) 
results in 


[«• V + a n (E)]<fi n . so{f) (x,n,E) = f”a mJ (E,E')<l> n ^(x,n,E')dE' 

+ f E (E, E ')<t> n . so(b) (x, -n, E')dE ' 


+ 


v r 


(30) 


® nk : iso(f ) 


+ e /; ® nk,iso(b ) (E,E')cf> kJor (x,-n,E')dE' 


for the forward component of the isotropic neutron flux. The backward moving isotropic neutron field satisfies a 
similar equation 


[-n-v + * n (E)]cfr n . S 0 {b) (x,n,E) = f~ a mJ {E,E'yt> n ^ h) (x,a,E ')dE' 

+ ?J>- Moif) (E,E')ct> kJor (x-n,E')dE' 

k 


T I* 

k 


J E a nk,iso(b) (E, E ')4> k for (x, Q,E')dE' 


( 31 ) 
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Note that (j) k j or (x, — f2, E ') = 0 in equations (30) and (31). Equations (30) and (31) represent a coupled set of 

differential equations that are solved efficiently using a Neumann series approach as shown by Slaba et al. [2010b]. 
This series is expected to converge quickly as o nn f(E,E') o nn ^(E,E } ). 

The isotropically produced neutrons initiate isotropic light ion production in collisions with the media 
resulting in 


B = J™a jn (E,E')<l> n . S0 (x,n,E')dE', (32) 

for which a semi-analytic solution can be efficiently implemented [Slaba et al. 2010b]. Equations (30) and (31) are 
solved simultaneously for no external fiuence incident on the media directed along Q (see Fig. 1). Note that the 
implicit assumption of equations (30) and (31) is that lateral diffusion from the ray 12 is fully compensated by 
diffusion from adjacent rays. This is strictly true for a flat plate slab with infinite lateral dimensions as will now be 
demonstrated and represents the use of periodic boundary conditions as an importance sampling method common to 
Monte Carlo simulation practice. 

These two methods (straight-ahead (N = 1) and bi-directional (N = 2)) are now compared to values 
simulated by the Monte Carlo codes Geant4, FFUKA, and PHITS for the Webber SPE spectrum normally incident 
on a 40 g/cm 2 slab of aluminum with infinite lateral dimensions. All Monte Carlo fiuence results presented herein 
are plotted as mean values with estimated statistical uncertainty shown as vertical error bars placed at the energy bin 
midpoint. The points are connected with lines to improve plot clarity. If not visible, vertical error bars are smaller 
than the symbols on the plot. Details of the Monte Carlo simulation setups are given in Appendix A. Results are 
shown in Fig. 8. The straight-ahead solution gives good agreement for charged components at most depths. Compare 
the left and right panes of Fig. 8. The neutron fiuence reaches a broad maximum across 20 to 35 g/cm 2 wherein the 
neutron fiuence is nearly isotropic at low energies and the straight-ahead (N = 1) fiuence approaches the bi- 
directional ( N = 2) fiuence from below before the maximum and over estimates beyond the maximum. On the back 
face, leakage is an important process leaving the straight-ahead fiuence as an over estimate at low neutron energies 
where the neutron fiuence is reduced by about half of its interior value due to leakage through the back surface 
(compare left and right panes of Fig. 8). The resonance structure seen in the Monte Carlo neutron fluences below 1 
MeV is a result of elastic cross section databases that have not been included in 3DHZETRN. 




10" 10" 1 10° IQ 1 10 2 10* 10‘ 3 H)' ] 10 c ' 1() ] HT IQ 3 

Kinetic energy (MeV) Kinetic energy (MeV) 

Fig. 8. Comparison of fiuence spectra at 35 g/cm 2 (left pane) and 40 g/cm 2 (right pane) in a 40 g/cm 2 aluminum slab 
(infinite in lateral dimensions) exposed to Webber SPE spectrum. 
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It is clear that the bi-directional solution for the semi-infinite slab where lateral leakage is negligible is a reasonable 
approximation. The success of the straight-ahead approximation in validation in human rated systems (Shuttle, ISS) 
clearly meets the requirements for a useful solution of the transport problem. Attention is now given to transport in 
finite objects in which lateral leakage plays an important role. 

3D Marching Procedures 

A sequence of approximations is now sought that will span the formalism from the simple bi-directional 
approximation to more complex propagation algorithms. With the demonstrated success of the bi-directional 
approximation, the solution is divided into forward and isotropic components, as before. The governing equation in 
the CSDA ignoring Coulomb scattering is given by equation (5). The particle flux is separated into forward and 
isotropic components, as in equation (20). Similarly, the double differential cross sections for neutron production are 
separated into forward and isotropic components given by equations (12) and (16). The charged particle production 
cross sections are assumed to take the form of equation (18). Upon substitution of the separated fluxes and 
production cross sections into equation (5), the forward components are obtained by solving 

B tt jJor (x,n,E)} = EJ e ^( E ’E%,fo r (x,n,E')dE', (33) 

k 

for charged particles (j ^ n ), and for neutrons (j = n ), 

[«• v + * n (E)]ct> nJor (x,n,E) = E Sy nkJor {E,E')ct> kJor ME')dE'. (34) 

k 

The space environment (external source) is nearly uniform in position and varies slowly in angle (if not 
isotropic), and for the forward flux, each directional component of the boundary condition may be propagated 
according to equations (33) and (34). Note that 12 enters equations (33) and (34) as a parameter and can be evaluated 
according to the boundary condition we have set as the incident flux from direction 12 0 ( see Fig 9). Once the forward 
flux is evaluated for direction 12 0 over the domain of x and E, one must yet solve for the induced isotropic neutron 
flux within the shield configuration. Coupling to the isotropically produced neutron field is taken as 


[n.v + a n (E)]cf> n!t jx,n,E) = J 00 / + 


where the isotropic neutron source is given by 


= Ej^n k , l J^ E '^0)Ajor(^ E ')dE' 

k 


(36) 


U-, 



Fig. 9. Geometry of 3D marching procedure. 
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Note that the isotropic neutron source is generated by collisions between the forward flux and target media, for 
which the forward flux is a function of only the penetration depth z (Fig. 9) along the direction 12 0 - This is denoted 
by writing the isotropic neutron source as 

x n , iso [<x),n,n 0 ,E} = z n .jx,n,n 0 ,E) , (37) 


where z(x) is the penetration depth along 12 0 to the point*. Equation (35) is now rewritten as 


[fi* V + a n (E)]cf> n . S o(x,n,E) = f 4 ^ m (E,E\n,nw n ^(x,n\EOdn'dE' + Xn , iso Hx),n,n 0 ,E] (38) 


Equation (38) also provides a basis for introducing anisotropic source terms as well. Note that while the forward 
propagating solution has its source fixed by the boundary condition, the isotropic solution (having no inbound 
fluence) is driven by internal sources generated by the collisions of the inbound fluence (as modified by the forward 
propagator in penetrating to * along direction 12 0 )- 

To solve equation (38), the neutron fluence <fi n iso (x, fl\E') must be known at every * from all directions IT 

approaching * from adjacent locations as was the case of equation (23). Although the use of approximation (25), 
based on the assumption that losses along the ray 12 are compensated by diffusion from adjacent rays, is a 
demonstrated useful assumption in a semi-infinite slab, it is an approximate overestimate in a sphere. Yet, using 
equation (25) in evaluation of equation (38) seems a reasonable approximation in most human rated vehicles leading 
to an extension of transport along an arbitrary ray 12 as a solution to the bi-directional equations. The forward 
component of the isotropic neutron field satisfies 


J »oo 

E °nnj ( E > E ')Kso(f) (x,fl,E')dE' 

+ iy n n^ E ')Kso^-^ E ') dE ' ( 39 ) 

+Xn,iso[ Z ( X )^’^0’E] , 

and the backward moving field satisfies 

[— fi» V + * n (E)]J) n . S0{b) (x,n,Ej = J^a nnJ (E,E')<P nttS0{b] (xME')dE' 

+ J^nnA E , E, Kiso(fM^ E ') dE ' ( 40 ) 

+x nMo Hx),-n,n 0 ,E] . 


It is clear that 12 enters equations (39) and (40) as a parameter and can be solved along an arbitrary number 
of rays used to construct a numerical representation of (/> n iso (&, 12, E) whose 12 dependence comes from distance to 

the boundary t( 12) and the penetration depth z along 12 0 . The isotropic neutron source of equation (38) effectively 
decouples the incident direction 12 0 from the neutron direction of propagation 12. 

These isotropically produced neutrons produce isotropic light ions in collisions with the media resulting in 
a light ion fluence given by 


B[4>j.Jx,n, E )] = (T jn (E, E ')4> niso {x, n,E')dE'. (41) 

Equations (39) and (40) are solved simultaneously for no external fluence incident on the media directed along 12 
(see Fig. 9). Nonetheless, efficient numerical methods already developed for the bi-directional model may be 
utilized. Note also that the implicit assumption of equations (39) and (40) is that lateral diffusion from the ray 12 is 
fully compensated by diffusion from adjacent rays. This has been demonstrated for slab geometry with results in 
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Fig. 8. The first order lateral leakage in this case is found by solving equations (39) and (40) along each £2/ over the 
ray limits {-£(-£2/), £(0/)} for a given set of N directions The next sections address convergence in the isotropic 
neutron and ion fields as a function of the number of directions ( N) chosen to represent the fields. Numerical 
solutions for equations (39) - (41) are the same as in the bi-directional solution and described by Slaba et al. 
[2010b]. The only difference here is that the isotropic source term in each case is computed as a function of 
penetration depth with respect to the incident source direction, £2 0 . 

The present transport formalism is related, but quite different from S N transport theory [Marchuk and 
Lebedev 1986]. The main difference between the two formalisms is that solutions along a set of directions 12/ are 
fully coupled to all other discrete directions in the S N treatment. While in the present formalism, coupling is only 
included through the isotropic source term and for transport along anti-parallel rays (i.e. 12/ and -12/). It is 
demonstrated below that the present formalism compares well to Monte Carlo simulations where such coupling is 
implicitly included. 


A Simple Geometry and Source 

Simple geometry and source orientation are chosen to allow efficient Monte Carlo simulation for 
comparison with the calculation of quantities with the above formalism. In this way, the solution methodology is 
verified and experience is gained in the formalism in preparation for connecting to the general geometric 
descriptions generated by the engineering design process for which Monte Carlo simulation is impractical. The 
geometry herein is taken as a simple sphere of aluminum of 40 g/cm 2 diameter for comparison with the appropriate 
solutions of the current formalism. Target points (detectors) are placed at depths along the z-axis (x = 0, y = 0). The 
external radiation environment (external source boundary condition) is assumed to be anti-parallel with the z-axis, 
uniform in the x-y plane, positioned above the sphere, and directed down onto the top of the sphere as shown in Fig. 
10 . 


^ ^ ^ Oo 



Fig. 10. Spherical geometry and external source orientation for benchmark comparisons (not all detectors shown). 


The object geometry is specified at a point x where field values are to be evaluated, about which the 
thickness of material t is known over an array of directions £2/ (i.e. t t = ^(£2/)). Note that each £2/ for i odd is backed 
by £2/+/ = -Q /. The remainder of the solution methodology lies in selection of the arbitrary directions of £2/ of 
propagation within the target material and evaluation of the penetration depth z for any transport depth x T along each 
£2/ for the given £2 0 at the boundary. The N directions are taken as the central ray of an equal division of the unit 
sphere such that A£2 = 4n / N. In this, the interest is in evaluating the omni-directional field function at an arbitrary 
location x about which the surrounding material is known as a thickness distribution of the various materials (in the 
present only homogenious aluminum). The evaluation will be made for a single direction of incidence £2 0 to 
highlight 3D transport effects, but this could include an array of incident directions for which the solutions would be 
integrated over all £2 0 values. 
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Fig. 11. Isotropic neutron source (E > 1 MeV) distribution induced by the Webber SPE spectrum in the spherical 

benchmark geometry. 


In that low-energy neutron transport is dominated by diffusive processes, it is anticipated that the solution 
will rapidly converge for low values of N (i.e. less than 20). The source of neutrons induced within the sphere by the 
Webber SPE is shown in Fig. 11. The resulting neutron fluence depends not only on the distance £(Q) to the 
boundary, but the intensity of the source along the ray in reaching the boundary. Clearly for SPEs, the most intense 
regions must be adequately represented in the domain of the solution. The next three sub-sections describe how the 
penetration depth, z, is determined for N= 1,2 and the general case (N> 2). 

N = 1, Straight-ahead approximation 

For each 12 0 , only one value of Q, is allowed and is taken to be that of the incident direction Q 0 - The 
transport distance x along 12 x and penetration depth z are the same, reducing complication in the formulation for N = 
1. In this case, leakage from the media is only through the distal face. The structure of this code is quite different 
than the original HZETRN code, although computational results are identical, and the original HZETRN code 
provides a verification test-case of the new 3DHZETRN code. For N = 1, the comparisons in the prior sections 
describing the verification and validation of HZETRN is fully applicable. If one restricts interest to N = 1 , then the 
HZETRN code is preferred in that computational efficiency is higher but the added overhead of 3DHZETRN for N 
= 1 is minimal. In the following sections, comparisons between 3DHZETRN and Monte Carlo simulations will also 
include results from the HZETRN code (N = 1) to highlight shortcomings of the straight-ahead approximation and 
improvements of the 3D corrections presented herein. 

N = 2, Bi-directional approximation 

The first step towards a full 3D code is the treatment of backward propagating particles. This was first 
accomplished by Clowdsley et al. [2000, 2001] and later improved by Slaba et al. [2010b] with demonstrated 
improved results [Slaba et al. 2010b, 2011b; Heinbockel et al. 2011a, 2011b]. In this case, propagation of the 
isotropic neutrons is represented in two directions by Qi = £2 0 and Q 2 = -£2 0 . Therefore, the penetration depth z is 
equal to x along fli and equal to ^(Oi)+ ^(Q 2 ) - * along Note that lateral leakage is not represented in this 
approximation although leakage at the forward and backward boundaries is so represented. The neutron transport in 
semi-infinite slab geometry has extensive Monte Carlo results for verification for which this two-stream 
approximation, ignoring lateral leakage, is most appropriate as demonstrated in Fig. 8. This approximation is 
appropriate for slab geometry and works well in applications in which the shield geometry is a collection of thin 
plates. The next step for treatment of lateral leakage is now discussed. 
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N> 2, N-stream approximation 

For general N, the transport formalism requires the penetration depth, z p , to be known as a function of 
transport depth, x T , along the propagation direction Q for each O 0 - In order to simplify algebraic manipulation and 
expressions, the center of the sphere is assumed to be positioned at the origin with the positive z-axis anti-parallel to 
the source direction, as shown in Fig. 10. Due to the simplified geometry and source orientation, an analytic 
expression for the penetration depth may be determined. The point x 0 may be written as 

x 0 = x — ft(t — x T ) , (42) 

where t is the known thickness along direction -Q from x to the boundary. The point x r may then be written in terms 
of x 0 as 


X T = X 0 - ct 0 z p . (43) 

The expression for z p is obtained by writing equation (43) in component form and substituting into the equation of a 
sphere with radius R = 20 g/cm 2 . This yields 


z p = x o' n o + p ? 2 - | x o | 2 + Ov^o ) 2 • ( 44 ) 

For preliminary testing and benchmark comparisons against Monte Carlo codes, the solutions for N = 6, 10, 
14, 18, and 22 are considered separately, in addition to the simple N = 1 and N = 2 cases discussed above. The 
distributions for all N are shown in Fig 12. Note that the N= 6, 14, and 22 cases include rays perpendicular to the z- 
axis (i.e. tangent to the sphere surface at the top and bottom detector locations). On the surface of the sphere, these 
rays correspond to zero shield thickness. This is corrected in the N= 10 and N = 18 distributions. The directional 
elements for each N are given in the Appendix B, along with a brief description of how the elements were 
determined. For each A, the omni-directional flux is obtained by integrating the fluxes along each ray direction after 
the transport procedure. 




Fig. 12. Directional components for N= 1, 2, 6, 10, 14, 18, and 22. 


Preliminary testing 

A first quantity of interest is the average path length at each detector location where the solution is 
evaluated to the sphere boundary. For a given N , the average path length is computed over the material thicknesses 
traversed from the target point to the boundary along all rays. This quantity is expected to converge rapidly at most 
locations with increasing N. As shown in Table 2, convergence is reasonably rapid for all locations except on the 
sphere surface. It is clear that the results are reasonably converged at N= 18, especially in comparison to N = 10, as 
will be further demonstrated for transport quantities. 
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Table 2. Average thickness traversed at different detector locations along z-axis of sphere. The left-most column is 
the distance from the top of the sphere to the detector. 


d (g/cm 2 ) 

Average thickness traversed (g/cm 2 ) 




N= 1 

N= 2 

N= 6 

N= 10 

N= 14 

N= 18 

N= 22 

0 

0.0 

20.0 

6.1 

10.4 

9.4 

10.1 

9.8 

0.3 

0.3 

20.0 

9.0 

10.9 

10.6 

10.7 

10.7 

5 

5.0 

20.0 

15.5 

15.6 

15.6 

15.7 

15.7 

10 

10.0 

20.0 

18.2 

18.2 

18.2 

18.2 

18.2 

15 

15.0 

20.0 

19.6 

19.6 

19.6 

19.6 

19.6 

20 

20.0 

20.0 

20.0 

20.0 

20.0 

20.0 

20.0 

25 

25.0 

20.0 

19.6 

19.6 

19.6 

19.6 

19.6 

30 

30.0 

20.0 

18.2 

18.2 

18.2 

18.2 

18.2 

35 

35.0 

20.0 

15.5 

15.6 

15.6 

15.7 

15.7 

40 

40.0 

20.0 

6.7 

10.4 

9.4 

10.1 

9.8 


The sequence of solutions for N = 1 to 22 is used as a test on the convergence of the solution method. In 
addition, a Monte Carlo benchmark is defined for each of the indicated methods. Benchmarks for N = 1 and N = 2 
are given in slab geometry in Fig. 8 to be later compared with N = 18 benchmarks in spherical geometry. 

The first quantity of interest is the dose associated with the Webber SPE as given in Table 3 for all N. Dose 
is mainly driven by the proton fluence that still has a large contribution from the penetrating primary protons that 
improves the convergence rate (especially at shallower depths near the top of the sphere). The dose equivalent, 
receiving a large contribution from the other light ions for which the quality factor is large, will require convergence 
of all the light ions to ensure an accurate value for dose equivalent. The values of dose equivalent for the Webber 
SPE are given in Table 4 for various numbers of rays in the 3D evaluations. Again, rapid convergence is observed 
with increasing N. As expected, the dose equivalent converges most slowly near the bottom of the sphere. A second 
look at the convergence of the particle fluence spectra are shown in Fig. 13. The convergence of the fluence spectra 
at 35 g/cm 2 is typical for an interior location for which N = 18 is clearly well converged. Even on the back surface 
of the sphere at 40 g/cm 2 the difference between N = 10 and N= 18 is very small in the graph even for alpha fluence; 
convergence is effectively obtained with N= 18. The rate of convergence of the alpha fluence mainly affected the 
convergence of the dose equivalent in Table 4. The N = 18 will be used in the Monte Carlo benchmarks discussed in 
the next section. 


Table 3. Dose (cGy/event) at different detector locations along z-axis of sphere for Webber SPE. The left-most 
column i s the distance from the top of the sphere to the detector. 


d (g/cm 2 ) 

Dose (cGy/event) 






N= 1 

N= 2 

N= 6 

N= 10 

N= 14 

N= 18 

N= 22 

0 

43703.70 

43703.77 

43703.73 

43703.84 

43703.79 

43703.82 

43703.80 

0.3 

1216.22 

1216.28 

1216.33 

1216.35 

1216.34 

1216.34 

1216.34 

5 

69.95 

69.98 

70.03 

70.02 

70.02 

70.02 

70.02 

10 

24.80 

24.81 

24.83 

24.82 

24.82 

24.82 

24.82 

15 

12.23 

12.23 

12.23 

12.23 

12.23 

12.23 

12.23 

20 

7.00 

6.99 

6.98 

6.98 

6.98 

6.98 

6.98 

25 

4.37 

4.36 

4.34 

4.34 

4.34 

4.34 

4.34 

30 

2.90 

2.88 

2.85 

2.85 

2.85 

2.85 

2.85 

35 

2.01 

1.98 

1.95 

1.95 

1.95 

1.95 

1.95 

40 

1.44 

1.41 

1.38 

1.38 

1.38 

1.38 

1.38 
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Table 4. Dose equivalent (cSv/event) at different detector locations along z-axis of sphere for Webber SPE. The left- 
mos t column is the distance from the top of the sphere to the detector. 


d (g/cm 2 ) 

Dose equivalent (cSv/event) 





N= 1 

N= 2 

N=6 

N = 10 

N = 14 

N = 18 

N= 22 

0 

737745.84 

737746.95 

737746.21 

737748.11 

737747.20 

737747.89 

737747.45 

0.3 

2871.17 

2872.15 

2873.07 

2873.31 

2873.10 

2873.17 

2873.12 

5 

109.95 

110.33 

111.15 

111.05 

111.02 

111.00 

111.00 

10 

37.94 

38.10 

38.45 

38.39 

38.36 

38.38 

38.38 

15 

19.01 

19.01 

19.04 

18.97 

19.00 

19.01 

18.99 

20 

11.33 

11.20 

10.97 

10.97 

10.98 

10.97 

10.98 

25 

7.54 

7.27 

6.85 

6.92 

6.90 

6.90 

6.91 

30 

5.42 

5.03 

4.58 

4.55 

4.59 

4.58 

4.58 

35 

4.15 

3.63 

3.18 

3.13 

3.13 

3.15 

3.15 

40 

3.34 

2.69 

2.26 

2.21 

2.19 

2.19 

2.21 




Fig. 13. Particle spectra at detector locations 35 g/cm 2 (left pane) and 40 g/cm 2 from top of aluminum sphere 
(diameter 40 g/cm 2 ) exposed to Webber SPE spectrum for different N in 3D transport solution. 


It is clear in the prior discussion that the main difference between N = 2 and N > 2 is the degree of lateral 
diffusion and the resulting leakage near boundaries. The ratio 


r N (E) 


1 </> n {E,N> 2) 

<t> n (E,N = 2) 


(45) 


is directly related to the representation of lateral losses for neutrons. The quantity r^{E) is fixed by the shape of the 
object, the scattering angular distribution, and the fidelity of the transport solution represented by N. The quantity 
r^{E) is shown in Fig. 14 for N = 6, 10, 14, 18, and 22. The value of zero at the highest energy is a result of the 
straight-ahead approximation for the forward component for which there is no lateral leakage. The largest relative 
difference between the N = 10 and N = 22 ratios is less than 3% at both depths, and the N = 18 solution appears 
converged based on the results in Tables 2-4 and Figs. 13 and 14. 
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Fig. 14. The neutron leakage ratio at detector locations 35 g/cm 2 (left pane) and 40 g/cm 2 (right pane) from top of 
aluminum sphere (diameter 40 g/cm 2 ) exposed to Webber SPE spectrum for different N in 3D transport solution. 
The largest relative difference between the N= 10 and N= 22 ratios is less than 3% at both depths. 


Monte Carlo Benchmarks 

In 2005, a systematic benchmark exercise was administered as a means of verifying various computational 
procedures resulting from the first phase of code development activity (1996-2004). Code verification uses 
established benchmarks as consisting of standard environments: the Webber SPE spectrum given in equation (19) 
and the 1977 solar minimum GCR environment. 

In earlier SPE benchmarking, the configuration was for a 30 cm semi-infinite slab of water shielded by 20 
g/cm 2 of either aluminum or iron utilizing only dose and dose equivalent (ICRP 26 quality factor [ICRP 1977]) as 
given in Figs. 5 and 6. For the moment, focus is mainly given to the results for the neutron spectra at the 
aluminum/water interface as shown in Fig. 15. It is seen from the figure that those codes with roots in the HETC 
code (HETC-HEDS and MCNPX) form a cluster of results including HZETRN at low energies ( E < 40 MeV). In 
distinction, Geant4 and FLUKA form a cluster below HZETRN. All codes show improved agreement above 100 
MeV. 



Fig. 15. Secondary neutron fluence at the interface between a 20 g/cm 2 aluminum shield and 30 g/cm 2 water target 

exposed to the Webber SPE (from Slaba et al. [2010b]). 
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It was generally concluded that HZETRN agreeds with the Monte Carlo codes as well as the various Monte 
Carlo codes agreed among themselves. There are currently no compelling reasons to attempt to use Monte Carlo 
methods in design and operational applications. However, it is further demonstrated below, the Monte Carlo codes 
provide an indispensable guide to development of the 3DHEZTRN code (that does provide a useful basis for space 
radiation protection). Focus is now shifted towards benchmarking 3DHZETRN, in which the added complications of 
finite geometry are addressed. Details associated with Monte Carlo simulations presented in Figs. 16-21 are given 
in Appendix A. One needs to be mindful of the intercode comparisons of Fig. 15 when viewing succeeding 
benchmarks. 

Webber SPE Benchmark 


The benchmark SPE spectrum is the Webber [1966] approximation of the February 23, 1956 SPE (not an 
accurate representation but a historically useful test spectrum, see Foelsche et al. [1974]) as given in equation (19). 
The induced neutron and proton fluence spectra were evaluated in a 40 g/cm 2 diameter aluminum sphere by Monte 
Carlo codes (Geant4, FEUKA, and PHITS) and 3DHZETRN with the results at 35 and 40 g/cm 2 depths shown in 
Fig 16. The proton fluences among the three codes are in reasonable agreement. The larger descrepancy lies with the 
neutron spectra. One interesting note is that 3DHZETRN agrees with Geant4 below 10 MeV and agrees with 
FEUKA above 10 MeV. FEUKA and 3DHZETRN exhibit a similar inflection in the neutron spectrum near 10 to 60 
MeV that may be related to the use of an early FEUKA approximation to the secondary neutron spectrum used in 
HZETRN and carried over into 3DHZETRN. Some of the differences in proton spectra on the back surface may be 
associated with the simplified assumptions of the straight-ahead approximation in low energy isotropic light ion 
transport. A principle difference in the three codes is the computational efficiency, as discussed in more detail later 
in this section. It is also worth noting that the 3DHZETRN (N= 18) solution is a substantial improvement over the 
HZETRN ( N = 1), indicating that lateral neutron leakage is more accurately described by the updated transport 
model. 




10 '' 10 ° 10 1 10 - 10 3 
Kinetic energy (MeV) 


Geant4 

FLUKA 

PIIITS 

3DHZETRN (N = 18) 
3DIIZETRN (N - 1) 


Fig. 16. Particle spectra at detector locations 35 g/cm 2 (left pane) and 40 g/cm 2 from top of aluminum sphere 
(diameter 40 g/cm 2 ) exposed to Webber SPE spectrum. 


The definition of the leakage factor given in equation (45) is specific to the 3DHZETRN formalism 
presented herein; however, an analogous quantity may also be computed for Monte Carlo results by comparing the 
ratio of the fluences in a slab (Fig. 8) to those obtained in a sphere (Fig. 16) provided the slab thickness and sphere 
diameters are identical. The leakage factor evaluated by three Monte Carlo codes (PHITS, FEUKA, and Geant4) are 
shown in comparison with 3DHZETRN (N= 18) at 35 g/cm 2 and 40 g/cm 2 in Fig. 17. The leakage is related to the 
angular dependence of the cross section. For example, the near zero leakage in the 3DHZETRN results above 200 
MeV is a result of the straight-ahead approximation. It is clear that the Monte Carlo results show significant 
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differences in leakage as determined mainly by differences in angular scattering cross sections. The 3DHZETRN 
results appear at a slight disadvantage with respect to Monte Carlo evaluations. It should also be clear that the 
improvement of the 3DHZETRN cross sections could result in more trustworthy results. This would not only 
improve leakage estimates but also give better representation of the resonance structure seen in the neutron fluences 
at lower energies. 




Gcant4 
FLUKA 
PH ITS 

3DHZETRN (N= 18) 


0.6 - 


10 “ 10 10 ' 

Kinetic energy (McV) 


10 - 10 3 10 '- 10 “ 10 10 1 

Kinetic energy (McV) 


Fig. 17. The neutron leakage ratio factor at detector locations 35 g/cm 2 (left pane) and 40 g/cm 2 (right pane) from 
top of aluminum sphere (diameter 40 g/cm 2 ) exposed to the Webber SPE spectrum. 
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The time required by the three codes to run these benchmarks is given in CPU seconds used and listed in 
the first column of Table 5. The 5 CPU seconds of 3DHZETRN allows for the development of affordable means of 
design optimization software on a commercial laptop while Geant4, FLUKA, and PHITS software requires a 
multiprocessing cluster environment to do a single design evaluation and would still be impractical for design 
optimization. The main limitation of the Monte Carlo methods is the number of particle histories required for 
reasonable convergence. It should be noted that source biasing, physics biasing, or other optimization strategies 
could significantly reduce the Monte Carlo runtimes shown in Table 5. However, the fully optimized run times 
would still be orders of magnitude larger than the 3DHZETRN run times and still require a multi-processor cluster 
computing environment. 


Table 5. Total CPU seconds required for benchmark calculations. Total number of histories ran in Monte Carlo 
s imulations are shown in parentheses. 



Webber SPE 

GCR hydrogen 

GCR helium 

GCR carbon 

GCR iron 

3DHZETRN (N= 1) 

2 

2 

2 

2 

2 

3DHZETRN (N= 18) 

5 

8 

8 

8 

8 

Geant4 

9xl0 7 (2xlO n ) 

1x10 s (2xlO n ) 

3xl0 8 (5xl0 9 ) 

2xl0 8 (2xl0 9 ) 

4xl0 8 (5xl0 8 ) 

FLUKA 

2xl0 8 (3xl0 n ) 

2xl0 8 (3xl0 n ) 

Ixl0 8 (lxl0 10 ) 

9xl0 7 (7xl0 9 ) 

2xl0 7 (5xl0 8 ) 

PHITS 

5xl0 10 (6xl0 8 ) 

7x10 9 (1x10 8 ) 

lxlO 9 (3xl0 8 ) 

2x10 8 (1x10 8 ) 

lxlO 8 (3xl0 8 ) 


Solar Minimum GCR benchmark. 

The present GCR benchmarks utilized the Badhwar-0 ’Neill 2010 [O’Neill 2010] model to represent the 1977 
solar minimum conditions (solar modulation parameter set to 475 MV). The hydrogen, helium, carbon, and iron 
spectra were considered separately as external source boundary conditions. The results for the neutron and proton 
fluence spectra within a 40 g/cm 2 sphere of aluminum for an incident solar minimum GCR hydrogen spectrum are 
shown in Fig. 18 at the two depths of 35 and 40 g/cm 2 . The GCR hydrogen results are qualitatively similar to results 
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obtained for the Webber SPE, only now the 3DHZETRN values favor Geant4 over the entire energy range except 
for the proton fluence below 100 MeV at the back surface (40 g/cm 2 ). The computation time used for transporting 
GCR hydrogen is slightly longer than for the Webber SPE transport in the same object. Still, the 3DHZETRN time 
is a very manageable 8 CPU seconds whereas the Monte Carlo codes take several orders of magnitude longer. 



Fig. 18. Particle spectra at detector locations 35 g/cm 2 (left pane) and 40 g/cm 2 from top of aluminum sphere 
(diameter 40 g/cm 2 ) exposed to the 1977 solar minimum GCR hydrogen spectrum. 
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Fig. 19. Particle spectra at detector locations 35 g/cm 2 (left pane) and 40 g/cm 2 from top of aluminum sphere 
(diameter 40 g/cm 2 ) exposed to the 1977 solar minimum GCR helium spectrum. The proton and alpha results have 

been scaled by 0.1 and 0.01 to improve plot clarity. 


The neutron, proton, and alpha spectra generated by GCR helium ions are shown in Fig. 19. A large 
contribution to the alpha fluence spectrum is from the incident helium ions. There is generally better agreement with 
Geant4 results than those generated by FLUKA. The computer resources used in these calculations are given in 
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Table 5. Similar results are shown for GCR carbon and iron ions in Figs. 20 and 21 with the corresponding resources 
used in Table 5. The light ion and neutron production nuclear database is more limited in HZETRN (and 
3DHZETRN) when the incident ion has Z > 2. First, the light ion and neutron fragments from the projectile are 
assumed to have the same velocity as the fragmenting ion itself. Second, the light ions and neutrons produced as 
target fragments have been ignored. Furthermore, the light ions are mostly treated within the straight-ahead 
approximation that mainly affects the light ion solution below 1 00 MeV/n. The spectral results will improve as these 
additional nuclear databases and corresponding transport procedures are added to the 3DHZETRN code. 
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Fig. 20. Particle spectra at detector locations 35 g/cm 2 (left pane) and 40 g/cm 2 from top of aluminum sphere 
(diameter 40 g/cm 2 ) exposed to the 1977 solar minimum GCR carbon spectrum. The proton and alpha results have 

been scaled by 0.1 and 0.01 to improve plot clarity. 
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Fig. 21. Particle spectra at detector locations 35 g/cm 2 (left pane) and 40 g/cm 2 from top of aluminum sphere 
(diameter 40 g/cm 2 ) exposed to the 1977 solar minimum GCR iron spectrum. The proton and alpha results have 

been scaled by 0.1 and 0.01 to improve plot clarity. 
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Although it is difficult to favor one code over another on the basis of the veracity of the solutions obtained, 
there are vast differences in computational time required to evaluate fluence spectra even in this simplest of 
geometries (homogeneous sphere of a single elemental material). Still, Monte Carlo methods provide an important 
test point for developing codes capable of meeting operational and especially design requirements. Computational 
speed has proven to be important in spaceflight validation in LEO where the time structure of the environment can 
be used to test various environmental components as was done using ISS and the Liulin instrument [Wilson et al. 
2007; Slaba et al. 2011a, 2013]. The main limitation on such studies remains the uncertainty in the environmental 
models (especially for the trapped environment) and uncertainty in nuclear cross sections. 

Conclusions 

The present version of 3DHZETRN shows improved agreement with Monte Carlo results in finite 
geometries where lateral leakage is non-negligible. Despite the improvements, the enhanced formalism and code has 
inherent limitations (aside from the simplified geometry of a homogeneous material) that carry over from the 
available bi-directional methods [Slaba et al. 2010b] in which the main share of light ions are still treated in the 
straight-ahead approximation as coupled to the forward propagating neutrons. Only the isotropically produced 
neutrons are first transported in 3D without the coupling to the isotropically produced light ion components. Once 
the isotropically produced neutrons are evaluated, the final coupling to the light ions is implemented (without further 
nuclear reactions) in a full 3D manner. Hence, only a fraction of the light ions have a 3D treatment, resulting in 
transport effects of unquantified magnitude. Future work will require extension to inhomogeneous media, complex 
geometries, and treatment of the cross-coupling of 3D neutron/light ion transport. Even so, the present expanded 
treatment of 3D effects shows great promise in the path to a fully 3D transport algorithm. 

The present benchmarks provide useful exercises for evaluation of computational procedures and 
atomic/nuclear database. It is now time for improving the physical description of the transport process that is still 
hampered by a lack of basic nuclear physics database developments. The codes for space design are still incomplete 
in both the physical description and corresponding implementation into computational procedures. Even so, the 
current development appears as state-of-the-art in computational procedures for spacecraft design. 
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Appendix A: Monte Carlo simulations 

In this section, a brief description is given of the Monte Carlo codes Geant4, FLUKA, and PHITS used to 
simulate particle spectra shown in Figs. 8 and 16 - 21. A discussion of the source, geometry, and simulation setup 
used for the codes is also given. 

Geant4 version 9.4.6 [Agostinelli et al. 2003] was used in the present calculations with the physics list 
QGSPBERTHP. This physics list utilizes the quark gluon string model (QGS) for high energy interactions with 
nucleons, pions, and nuclei. Post-interaction nuclear de-excitation is handled by the precompound (P) model. The 
Bertini cascade (BERT) model is used for interactions below 10 GeV. A high-precision (HP) tracking model is used 
for neutrons below 20 MeV. Nucleus-nucleus collisions are represented by the native Quantum Molecular Dynamics 
(QMD) model [Koi 2008]. More information about the chosen physics list can be found at the Geant4 website 
[Geant4 2012a]. The QGSP BERT HP physics list has been recommended by the Geant4 collaboration for high 
energy physics applications [Geant4 2012b] and has been used, with some slight variation, in other space radiation 
related studies [Bemabeu and Casanova 2007, Hayatsu et al. 2008, Martinez and Kingston 2012, Slaba et al. 2013]. 

FLUKA is a Monte Carlo program that contains fully integrated physics and performs the transport of 
elementary particles and ions through materials. It has the ability to transport and calculate interactions with all 
elementary hadrons, light and heavy ions, and electrons and photons over an energy range which extends up to 10 4 
TeV for all particles, and down to thermal energies for neutrons [Battistoni et al. 2007]. The nuclear models are 
integrated into the software for all particles and energies, with the exception of neutrons below 20 MeV, where 
standard international evaluated data files are used. This code has been extensively benchmarked against available 
accelerator and cosmic ray experimental data, at beam energies as low as a few MeV and as large as cosmic ray 
energies [Aarnio et al. 1993; Andersen et al. 2004]. 
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PHITS has been used extensively in space radiation transport. It simulates ionization through transport 
processes and uses the mean free path to determine when collisions occur [Niita et al. 2006, Sihver et al. 2007]. 
Ionization transport can include angle and energy straggling, showing good agreement for Bragg peaks and 
fragmentation tails [Niita et al. 2006]. Tabulated nuclear cross section data are used for neutron interactions below 
20 MeV, and evaluated data are used for gamma and electron transport below 1 GeV [Niita et al. 2006]. High 
energy neutron and heavy charged particle interactions are simulated through the Jet AA Microscopic transport 
model (JAM), which addresses hadron-hadron interactions, and the Japan Atomic Energy Research Institute 
(JAERI) Quantum Molecular Dynamics model (JQMD), which addresses hadron-nucleus and nucleus-nucleus 
interactions, up to energies of 200 GeV/n [Niita et al. 2006, Sato et al. 2006, Sihver et al. 2007]. Both JAM and 
JQMD are used for the dynamic portion of the nuclear reaction, and once the dynamics are completed, the General 
Evaporation Model is used to address nuclear de-excitation through particle emission in a statistical manner [Niita et 
al. 2006, Sihver et al. 2007]. The total nucleus-nucleus cross sections used in PHITS are the same as those that have 
been used in versions of HZETRN [Niita et al. 2006]. 

In Fig. 8, the geometry setup was for a 40 g/cm 2 aluminum slab with infinite lateral dimensions and the 
external source was applied normally to the front face of the slab. Particle spectra were tallied at a thin (0.027 g/cm 2 ) 
aluminum detector layer placed 35 g/cm 2 and 40 g/cm 2 from the front of the slab. In Figs. 16 - 21, the geometry 
setup was for a sphere with diameter 40 g/cm 2 , and the external source was applied uniformly onto the top of the 
sphere (covering the entire sphere) parallel with the z-axis as shown in Fig. 10. Particle spectra were tallied at 
spherical aluminum detectors (radius 0.675 g/cm 2 ) placed along the z-axis down through the sphere. Convergence 
tests on the detector radius were performed to arrive at the final choice. In general, the chosen detector radius is 
large enough to allow reasonable tally statistics, but small enough to effectively prevent fragmentation events from 
occurring within the detector. 

For all of the GCR simulations, the source energy spectrum was sampled directly from the input spectrum 
(i.e. no source energy biasing). For the SPE simulations, the source energy spectrum was biased so that source 
sampling was performed from a uniform distribution covering 0.01 MeV to 2500 MeV. This biasing technique 
places greater emphasis on the higher energy portion of the SPE spectrum improving tally statistics for secondary 
neutrons. 

Appendix B: Evaluation of Q, for N= 6, 10, 14, 18, 22 

The Q, are herein evaluated with £2 0 as the top hemispheric pole so that = Q 0 , while the bottom pole is at 
Q 2 = -^o and are the central rays of the two polar regions. The remaining directions denoted by Q, for i = 3...N are 
symmetrically distributed about the mid-latitudes. Generally in application, the Q,’s must be rotated to the body 
coordinates of the transport media. This is accomplished by a rotation matrix R(6^, </>o) determined from the polar 
vector Q 0 relative to the body coordinate system (note, 6 denotes co-latitude with values on 0 to 7i, and ^ denotes 
longitude running from 0 to 2n). The regional boundaries are chosen such that 

f dfi =47 T / N, (46) 

where AQ* denotes the z th region, N is the number of regions and is related to the number of Nq segments of Acos 6 t 
and the number of N# segments Afc plus the two polar regions such that 

N = N^N e + 2, (47) 

and defined such that 

F Af2 » = F A cos e i ■ ( 48 ) 

i = 1 i = 1 

The polar caps are bound by 


cos 0 1 = 1 — 2 / N , 


(49) 
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cos# 2 = — 1 + 2 / N , 


(50) 


and the corresponding central rays Qi and Q 2 lie along the ±z axis. The middle latitudes are divided into N^o 
regions extending from the polar cap to the equator in each hemisphere. The corresponding longitudinal boundaries 
are set by the N$ such that the 

A^=27t /N+. (51) 

Note that N# is generally greater than or equal to 4. The remaining boundaries in latitude are given as 

A cos 0 i = (cos 9 1 — cos # 2 ) / N 0 . (52) 

For N = 6, there is only one midpoint in latitude taken as 

cos#- = (cos0 1 + cos# 2 ) / 2 = 0, (53) 

for i> 2. The four longitudinal points are given as 

(f) = {0,7r / 2,7t, 3tt / 2}. (54) 

The directional elements for each solid angle region are given by 

ft !• = { cos (/)• sin , sin 0. sin 6 ] , cos 6 l } , (55) 

and are given for Af = 6 in Table B.l. Similar results for other N used in this paper are given in Tables B.2 - B.5. 

Table B.l Directional e lements of the Q f - for N= 6. 


i 

1 

2 

3 

4 

5 

6 

e x 

0 

0 

1 

-1 

0 

0 

e v 

0 

0 

0 

0 

1 

-1 

e z 

1 

-1 

0 

0 

0 

0 


Table B.2 Directional 

elements of the Q, for N = 

10. 






i 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 


0 

0 

0.9165 

-0.9165 

0 

0 

-0.9165 

0.9165 

0 

0 

e y 

0 

0 

0 

0 

0.9165 

-0.9165 

0 

0 

-0.9165 

0.9165 

e z 

1 

-1 

0.4 

-0.4 

0.4 

-0.4 

0.4 

-0.4 

0.4 

-0.4 

Table B.3 Directional elements of the Q* for N = 

14. 







i 

1 

3 

5 

7 

9 

11 

13 


0 

0.8207 

0 

-0.8207 

0 

1 

0 

e v 

0 

0 

0.8207 

0 

-0.8207 

0 

1 

e z 

1 

0.5714 

0.5714 

0.5714 

0.5714 

0 

0 

i 

2 

4 

6 

8 

10 

12 

14 

&x 

0 

-0.8207 

0 

0.8207 

0 

-1 

0 

e v 

0 

0 

-0.8207 

0 

0.8207 

0 

-1 

e z 

-1 

-0.5714 

-0.5714 

-0.5714 

-0.5714 

0 

0 
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Table B.4 Directional elements of the Q/ for N = 18. 


i 

1 

3 

5 

7 

9 

11 

13 

15 

17 

e x 

0 

0.7454 

0 

-0.7454 

0 

0.9750 

0 

-0.9750 

0 

e v 

0 

0 

0.7454 

0 

-0.7454 

0 

0.9750 

0 

-0.9750 


1 

0.6667 

0.6667 

0.6667 

0.6667 

0.2222 

0.2222 

0.2222 

0.2222 

i 

2 

4 

6 

8 

10 

12 

14 

16 

18 


0 

-0.7454 

0 

0.7454 

0 

-0.9750 

0 

0.9750 

0 

e v 

0 

0 

-0.7454 

0 

0.7454 

0 

-0.9750 

0 

0.9750 

e z 

-1 

-0.6667 

-0.6667 

-0.6667 

-0.6667 

-0.2222 

-0.2222 

-0.2222 

-0.2222 


Table B.5 Directional elements of the Q / for N = 22. 

i 1 

3 

5 

7 

9 

11 

13 

15 

17 

19 

21 

e x 0 

0.6863 

0 

-0.6863 

0 

0.9315 

0 

-0.9315 

0 

1 

0 

e v 0 

0 

0.6863 

0 

-0.6863 

0 

0.9315 

0 

-0.9315 

0 

1 

e z 1 

0.7273 

0.7273 

0.7273 

0.7273 

0.3636 

0.3636 

0.3636 

0.3636 

0 

0 

i 2 

4 

6 

8 

10 

12 

14 

16 

18 

20 

22 

e x 0 

-0.6863 

0 

0.6863 

0 

-0.9315 

0 

0.9315 

0 

-1 

0 

e y 0 

0 

-0.6863 

0 

0.6863 

0 

-0.9315 

0 

0.9315 

0 

-1 

e * - 1 

-0.7273 

-0.7273 

-0.7273 

-0.7273 

-0.3636 

-0.3636 

-0.3636 

-0.3636 

0 

0 


References 

Aamio, P.A., Fasso, A. Ferrari, A., Moehring, H.J., Ranft, J., Sala, P.R., Stevenson, G.R., Zazula, J.M., FLUKA: 
Hadronic benchmarks and applications. Proceedings of the MC93 International Conference on Monte Carlo 
Simulation in High Energy and Nuclear Physics, 1993. 

Adamczyk, A.M., Norman, R.B., Sriprisan, S.I., Townsend, L.W., Norbury, J.W., Blattnig, S.R., Slaba, T.C., 
NUCFRG3: Light ion improvements to the nuclear fragmentation model. Nucl. Instrum. & Methods A 678: 21-32; 
2012. 

Agostinelli, S. et al., Geant4— a simulation toolkit. Nucl Instrum. & Methods A 506: 250-303; 2003. 

Alsmiller, R.G., Irving, D.C., Kinney, W.E., Moran, H.S., The validity of the straightahead approximation in space 
vehicle shielding shielding studies. In Second Symposium on Protection Against Radiations in Space , Arthur Reetz, 
ed. NASA SP-71, pp. 177-181, 1965. 

Andersen, V., Ballarini, F., Battistoni, G., Campanella, M., Carboni, M., Cerutti, F., Empl, A., Fasso, A., Ferrari, A., 
Gadioli, E., Garzelli, M.V., Lee, K., Ottolenghi, A., Pelliccioni, M., Pinsky, L.S., Ranft, J., Roesler, S., Sala, P.R., 
Wilson, T.L., The FLUKA code for space applications: recent developments. Adv. Space Res. 34: 1338-1346; 2004. 

Armstrong, T.W., Colburn, B.L., Predictions of secondary neutrons and their importance to radiation effects inside 
the International Space Station. Radiat. Meas. 33: 229-234; 2001. 

Badhwar, G.D., Patel, J.U., Cuccinotta, F.A., Wilson, J.W., Measurements of the secondary particle energy spectra 
in the space shuttle, Rad. Meas., 24: 129-138; 1995. 

Badhwar, G.D., O’Neill, P.M., Cucinotta, F.A., In-flight radiation measurements on STS-60. Radiat. Meas. 26:17- 
34; 1996. 

Badhwar, G.D., Atwell, W., Badavi, F.F., Yang, T.C., Space radiation absorbed dose distribution in human 
phantom. Radiat. Meas. 33: 235-241; 2001. 


27 



Battistoni, G., Muraro, S., Sala, P.R., Cerutti, F., Ferrari, A., Roesler, S.,Fasso, A., Ranft, J., The FLUKA code: 
Description and benchmarking. Proceedings of the Hadronic Shower Simulation Workshop 2006 , 896 : 31-49; 2007. 

Bemabeu, J. and Casanova, I., Geant4-based radiation hazard assessment for human exploration missions. Adv. 
Space Res. 40 : 1368-1380; 2007. 

Clowdsley, M.S., Heinbockel, J.H., Kaneko, FL, Wilson, J.W., Singleterry, R.C., Shinn, J.L., A comparison of the 
multigroup and collocation methods for solving the low-energy neutron Boltzmann equation. Can. J. Phys. 78 : 45- 
56; 2000. 

Clowdsley, M.S., Wilson, J.W., Shinn, J.L., Badavi, F.F., Heinbockel, J.H., Atwell, W., Neutron environment 
calculations for low Earth orbit. SAE 01ICES2327, 2001. 

Cucinotta, F.A. and Dubey, R.R., Alpha-cluster description of excitation energies in 12 C( 12 C,3a)X at 2.1 A GeV, 
Phys. Rev. C 50 : 979-984; 1994. 

Cucinotta, F.A., Wilson, J.W., Shinn, J.L., Tripathi, R.K., Assessment and requirements of nuclear reaction data 
bases for GCR transport in the atmosphere and structures, Adv. Space Res., 21 ( 12 ): 1753-1762; 1998. 

Cucinotta, F.A., Wilson, J.W., Saganti, P., Hu, X. Kim, M.Y., Cleghorn, T., Zeitlin, C.E., Tripathi, R.K., Isotopic 
dependence of GCR fluence behind shielding. Radiat. Meas. 41 ( 9 ): 1235-1249; 2007. 

Fasso, A., Ferrari, A., Ranft, J., Sala, P.R., FLUKA: A multi-particle transport code, CERN-2005-10, INFN/TC 
05/11, SLAC-R-773, 2005. 

Foelsche, T., Mendell, R.B., Wilson, J.W., Adams, R.R., Measured and calculated neutron spectra and dose 
equivalent rates at high altitudes: Relavance to SST operations and space research. NASA TN D-7715, 1974. 

Geant4 collaboration. Geant4 Reference Physics Lists. Website: 

http://www.geant4.org/geant4/support/physicsLists/referencePL/referencePL.shtml (20 1 2a). 

Geant4 collaboration. Geant4 Recommended Physics Lists. Website: 
http://www.geant4.org/geant4/support/physicsLists/referencePL/useCases.shtml (20 1 2b). 

Golovchenko, A.N., Skvarc, J., Ilic, R., Sihver, L., Bamblevski, V.P., Tretyakova, S.P., Schardt, D., Tripathi, R.K., 
Wilson, J.W., R. Bimbot, R., Fragmentation of 200 and 244 MeV/u carbon beams in thick tissue-like absorbers. 
Nucl. Instrum. & Methods B, 159 : 233-240; 1999. 

Hayatsu, K., Hareyama, M., Kobayashi, S., Yamashita, N., Miyajima, M., Sakurai, K., Hasebe, N., Radiation doses 
for humans exposed to galactic cosmic rays and their secondary products on the lunar surface. Bio. Sci. in Space 22 : 
59-66; 2008. 

Heinbockel, J.H., Feldman, G.A., Wilson, J.W., Singleterry, R.C., Leakeas, C.L., Clowdsley, M.S., Solutions to the 
low energy neutron Boltzmann equation for space applications. SAE 03ICES339, 2003. 

Heinbockel, J.H., Slaba, T.C., Blattnig, S.R., Tripathi, R.K., Townsend, L.W., Handler, T., Gabriel, T.A., Pinsky, 
L.S., Reddell, B., Clowdsley, M.S., Singleterry, R.C., Norbury, J.W., Badavi, F.F., Aghara, S.K., Comparison of the 
transport codes HZETRN, HETC and FLUKA for a solar particle event. Adv. Space Res. 47 : 1079-1088; 2011a. 

Heinbockel, J.H., Slaba, T.C., Tripathi, R.K., Blattnig, S.R., Norbury, J.W., Badavi, F.F., Townsend, L.W., Handler, 
T., Gabriel, T.A., Pinsky, L.S., Reddell, B., Aumann, A.R., Comparison of the transport codes HZETRN, HETC and 
FLUKA for galactic cosmic rays. Adv. Space Res. 47:1089-1 105; 2011b. 

Hugger, C.P., Qualls, G.D., Wilson, J.W., Cucinotta, F.A., Shavers, M.R., Zapp, N., ISS as a platform for 
environmental evaluation. AIP Conf. Proc. 654: 1011-1017; 2003. 


28 



ICRP, Recommendations of the ICRP. ICRP Publication 26, 1977. 

Koi, T., New native QMD code in Geant4. IEEE Nuclear Science Symposium Conference Record , 2008. 

Lamkin, S.L., A theory for high-energy nucleon transport in one-dimension. Masters Thesis, Old Dominion 
University, June 1974. 

Nealy, J.E., Cucinotta, F.A., Wilson, J.W., Badavi, F.F., Dachev, Ts.P., Tomov, B.T., Walker, S.A., de Angelis, G., 
Blattnig, S.R., Atwell, W., Pre-engineering spaceflight validation of environmental models and the 2005 HZETRN 
simulation code. Adv. Space Res. 40 ( 11 ): 1592-1610; 2006. 

Martinez, F.M. and Kingston, J., Space radiation analysis: Radiation effects and particle interaction outside the 
earth’s magnetosphere using GRAS and GEANT4. Acta Astronautica 72 : 156-164; 2012. 

Niita, K., Sato, T., Iwase, H., Nose, H., Nakashima, H., Sihver, F., 2006 PHITS - A particle and heavy ion transport 
code system, Rad. Meas ., 41 : 1080 - 1090; 2006. 

Marchuk, G.I. and Febedev, V.I., Numerical Methods in the Theory of Neutron Transport , Harwood Academic 
Publishers, 1986. 

Norbury, J.W., Adamczyk, A., Differential cross sections for electromagnetic dissociation. Nucl. Instrum. & 

Methods B 254 : 177-181; 2007. 

Norman, R.B., Blattnig, S.R., De Angelis, G., Badavi, F.F., Norbury, J.W., Deterministic pion and muon transport in 
Earth’s atmosphere. A dv. Space Res. 50 : 146-155; 2012. 

Norman, R.B., Slaba, T.C., Blattnig, S.R., An extension of HZETRN for cosmic ray initiated electromagnetic 
cascades. Adv. Space Res. 51 : 2251-2260; 2013. 

O’Neill, P.M., Badhwar-O’Neill galactic cosmic ray flux model - Revised. IEEE Trans. Nuc. Sci., 57 : 3148-3153; 

2010. 

Pinsky, F., Carminati, F., Ferrari, A., Simulation of space shuttle neutron measurements with FFUKA. Radiat. 

Meas. 33 : 335-339; 2001. 

Ranft, J., The FFUKA and KASPRO hadronic cascade codes. Computer Techniques in Radiation Transport and 
Dosimetry, W. R. Nelson and T.M. Jenkins, eds, Plenum Press, pp. 339-371, 1980 

Sato, T., Niita, K., Iwase, H., Nakashima, H., Yamaguchi, Y., Sihver, F., 2006 Applicability of particle and heavy 
ion transport code PHITS to the shielding design of spacecrafts, Rad. Meas., 41 : 1 142 - 1 146; 2006. 

Sato, T., Niita, K., Matsuda, N., Hashimoto, S., Iwamoto, Y., Noda, S., Ogawa, T., Iwase, H., Nakashima, H., 
Fukahori, T., Okumura, K., Kai, T., Chiba, S., Furuta T., Sihver, F., Particle and heavy ion transport code system 
PHITS, version 2.52, J. Nucl. Sci. Technol. 50 : 913-923; 2013. 

Schimmerling, W., Rapkin, M., Wong, M., Howard, J., The propogation of relativistic heavy ions in multielement 
beam lines. Medical Phys. 13 : 217-228; 1986. 

Scott, W.W., Alsmiller, R.G., Comparisons of results obtained with several proton penetration codes. ORNF-RSIC- 
17. 1967. 

Shavers, M.R., Frankel, K., Miller, J., Schimmerling, W., Townsend, F.W., Wilson, J.W., The fragmentation of 
670 A MeV Neon-20 as a function of depth in water. III. Analytic Multigeneration Transport Theory. Radiat. Res., 
134 ( 1 ): 1-14; 1993. 


29 



Shinn, J.L., Cucinotta, F.A., Simonsen, L.C., Wilson, J.W., Badavi, F.F., Badhwar, G.D., Miller, J., Zeitlin, C., 
Heilbronn, L., Tripathi, R.K., Clowdsley, M.S., Heinbockel, J.H., Xapsos, M.A., Validation of a comprehensive 
space radiation transport code, IEEE Trans. Nuc. Sci ., 45 : 2711-2719; 1998. 

Sihver, L., Mancusi, D., Sato, T., Niita, K., Iwase, H., Iwamoto, Y., Matsuda, N., Nakashima, H., Sakamoto, Y., 
2007 Recent developments and benchmarking of the PHITS Code. Adv. Space Res., 40 : 1320 - 1331; 2007. 

Singleterry, R.C., Blattnig, S.R., Clowdsley, M.S., Qualls, G.D., Sandridge, C.A., Simonsen, L.C., Slaba, T.C., 
Walker, S.A., Badavi, F.F., Spangler, J.L., Aumann, A.R., Zapp, E.N., Rutledge, R.D., Lee, K.T., Norman, R.B., 
Norbury, J.W., OLTARIS: On-line tool for the assessment of radiation in space. Acta Astronautica 68 : 1086-1097; 
2011. 

Slaba, T.C., Blattnig, S.R., Badavi, F.F., Faster and more accurate transport procedures for HZETRN. J. Comp. 
Phys. 229 : 9397-9417; 2010a. 

Slaba, T.C., Blattnig, S.R., Aghara, S.K., Townsend, L.W., Handler, T., Gabriel, T.A., Pinsky, L.S., Reddell, B., 
Coupled neutron transport for HZETRN. Radiat. Meas. 45 : 173-182; 2010b. 

Slaba, T.C., Blattnig, S.R., Badavi, F.F., Stoffle, N.N., Rutledge, R.D., Lee, K.T., Zapp, E.N., Dachev, Ts. P., 
Tomov, B.T., Statistical validation of HZETRN as a function of vertical cutoff rigidity using ISS measurements. 
Adv. Space Res. 47 : 600-610; 2011a. 

Slaba, T.C., Blattnig, S.R., Clowdsley, M.S., Variations in lunar neutron dose estimates. Rad. Res. 176 : 827-841; 
2011b. 

Slaba, T.C., Blattnig, S.R., Reddell, B., Bahadori, A., Norman, R.B., Badavi, F.F., Pion and electromagnetic 
contribution to dose: Comparisons of HZETRN to Monte Carlo results and ISS data. Adv. Space Res. 52 : 62-78; 
2013. 

Tai, H., Bichsel, H., Wilson, J.W., Shinn, J.L., Cucinotta, F.A., Badavi, F.F., Comparison of stopping power and 
range data bases for radiation transport study, NASA TP-3644, 1997. 

Tripathi, R.K., Townsend, L.W., Kahn, F., Role of intrinsic width in fragment momentum distributions in heavy ion 
collisions, Phys. Rev. C. 48(4): R1775-R1777; 1994. 

Tweed, J., Walker, S.A., Wilson, J.W., Cucinotta, F.A., Tripathi, R.K., Blattnig, S., Mertens, C.J., Computational 
methods for the HZETRN code. Adv. Space Res. 35: 194-201; 2006a. 

Tweed, J., Walker, S.A., Wilson, J.W., Tripathi, R.K., Cucinotta, F.A., Badavi, F.F., An improved Green’s function 
code for HZE transport. SAE/ICES paper 2006-01-2147, 2006b. 

Walker, S.A., Tweed, J., Wilson, J.W., Cucinotta, F.A., Tripathi, R.K., Blattnig, S., Zeitlin, C., Heilbronn, L., 
Miller, J., Validation of the HZETRN code for laboratory exposures with 1 A GeV iron ions in several targets. Adv. 
Space Res. 35: 202-207; 2006. 

Webber, W.R., An evaluation of solar-cosmic-ray events during solar minimum. D2-84274-1, Boeing Co. 1966. 

Wilson, J.W., Khandelwal, G.S., Proton dose approximation in convex geometry. Nucl. Tech. 23: 298-305; 1974. 

Wilson, J.W. and Lamkin, S.L, Perturbation theory for charged-particle transport in one dimension. Nucl. Sci. & 
Eng. 57(4): 292-299; 1975. 

Wilson, J.W., Analysis of the theory of high-energy ion transport. NASA TN D-8381, 1977. 

Wilson, J.W., Townsend, L.W., Bidasaria, H.B., Schimmerling, W., Wong, M. Howard, J., 20-Ne depth-dose 
relations in water. Health Phys. 46 : 1101-1111; 1984. 


30 



Wilson, J.W. and Badavi, F.F., Methods of galactic heavy ion transport. Radiat. Res. 108 : 231-237; 1986. 

Wilson, J.W., Townsend, L.W., Badavi, F.F., A semiempirical nuclear fragmentation model. Nucl. Instrum. & 
Methods B 18 : 225-231; 1987a. 

Wilson, J.W., Townsend, L.W., Badavi, F.F., Galactic cosmic ray propagation in earth’s atmosphere. Radiat. Res. 
109 : 173-183; 1987b. 

Wilson, J.W., Chun, S. Y., Buck, W.W., Townsend, L.W., High energy nucleon data bases. Health Phys ., 55: 817- 
819; 1988. 

Wilson, J.W., Townsend, L.W., Schimmerling, W., Khandelwal, G.S., Khan, F., Nealy, J.E., Cucinotta, F.A., 
Simonsen, L.C., Shinn, J.L., Norbury, J.W., Transport methods and interactions for space radiations, NASA RP- 
1257,1991. 

Wilson, J. W., Shavers, M.R., Badavi, F.F., Miller, J., Shinn, J.L., Costen, R.C., Nonperturbative methods in hze 
propagazition. Radiat. Res. 140 : 241-244; 1994a. 

Wilson, J.W., Townsend, L.W., Shinn, J.L., Badavi, F.F., Lamkin, S.L., Galactic cosmic ray transport methods: 
Past, present, and future. Adv. Space Res. 14 : 841-852, 1994b. 

Wilson, J.W., Shinn, J.L., Townsend, L.W., Tripathi, R.K., Badavi, F.F., Chun, S.Y., NUCFRG2: A semiempirical 
nuclear fragmentation model. Nucl. Instrum. & Methods B 94 : 95-102; 1994c. 

Wilson, J.W ., Tripathi, R. K., Cucinotta , F. A., Shinn, J. L., Badavi, F. F., Chun, S. Y., Norbury, J. W., Zeitlin, 

C. J., Hielbronn, L. H., Miller, J., NUCFRG2, an evaluation of the semiempirical nuclear fragmentation database, 
NASA TP-3533, 1995. 

Wilson, J.W, Miller, J.M., Konradi, A., Cucinotta, F.A., Shielding strategies for human space exploration. NASA 
CP-3360, 1997. 

Wilson, J.W., Cucinotta, F.A., Tai, H., Shinn, J.L., Chun, S.Y., Tripathi, R.K., Sihver, L., Transport of light ions in 
matter. Adv. Space Res. 21(12): 1763-1771; 1998. 

Wilson, J.W., Tweed, J., Heinbockel, J.H., Tripathi, R.K., Cucinotta, F.A., Mertens, C.J., Advances in space 
radiation transport codes. 14 th Space Radiation Health Investigators’ Workshop, League City, TX; 2003a. 

Wilson, J.W., Nealy, J.E., De Angelis, G., Badavi, F.F., Hugger, C.P., Cucinotta, F.A., Kim, M.Y., 
Dynamic/anisotropic low earth orbit environmental models. AIAA Space Conference paper 2003-6221, 2003b. 

Wilson, J.W., Tripathi, R.K., Qualls, G.D., Cucinotta, F.A., Prael, R.E., Norbury, J.W., Heinbockel, J.H., Tweed, J., 
Space radiation transport methods development, Adv. Space Res. 34 ( 6 ): 1319-27; 2004a. 

Wilson, J.W., Cucinotta, F.C., Schimmerling, W.S., Emerging radiation health-risk mitigation technologies. AIP 
Conf. Proc. 699 , 913-924; 2004b. 

Wilson, J.W., Tripathi, R.K., Mertens, C.J., Blattnig, S.R., Clowdsley, M.S., Verification and validation: High 
charge and energy (HZE) transport codes and future development. NASA TP-2005-213784, 2005. 

Wilson, J.W., Tripathi, R.K., Badavi, F.F., Cucinotta, F.A., Standardized radiation shield design method: 2005 
HZETRN, SAE/ICES paper 2006-01-2109, 2006. 

Wilson, J.W., Nealy, J.E., Dachev, Ts.P., Tomov, B.T., Cucinotta, F.A., Badavi, F.F., De Angelis, G., Atwell, W., 
Leutke, N., Time serial analysis of the induced LEO environment within the ISS 6A. Adv. Space Res. 40 ( 1 1 ): 1562- 
1570; 2007. 


31 



Zeitlin, C., Heilbronn, F., Miller, J., Rademacher, S.E., Borak, T., Carter, T.R., Frankel, K.A., Schimmerling, W., 
Stronach, C.E., Heavy fragment production cross sections from 1.05 GeV/nucleon 56Fe in C, Al, Cu, Pb, and CH 2 
targets. Phys. Rev. C 56 : 388-397; 1997. 

Zeitlin, C., Miller, J., Guetersloh, S., Heilbronn, L., Fukumura, A., Iwata, Y., Murakami, T., Blattnig, S., Norman, 
R.. Mashnik, S., Fragmentation of 14 N, 16 0, 20 Ne, and 24 Mg nuclei at 290 to 1000 MeV/Nucleon. Phys. Rev. C 83 : 
034909; 2011. 


32 



REPORT DOCUMENTATION PAGE 


Form Approved 
OMB No. 0704-0188 


l 

The public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, 
gathering and maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this 
collection of information, including suggestions for reducing this burden, to Department of Defense, Washington Headquarters Services, Directorate for Information Operations and 
Reports (0704-0188), 1215 Jefferson Davis Highway, Suite 1204, Arlington, VA 22202-4302. Respondents should be aware that notwithstanding any other provision of law, no person 
shall be subject to any penalty for failing to comply with a collection of information if it does not display a currently valid OMB control number. 

PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS. 


1. REPORT DATE (DD-MM-YYYY) 2. REPORT TYPE 3. DATES COVERED (From - To) 

01-05-2014 T echnical Publication 

4. TITLE AND SUBTITLE 

A 3DHZETRN Code in a Spherical Uniform Sphere with Monte Carlo 
Verification 

5a. CONTRACT NUMBER 

5b. GRANT NUMBER 

5c. PROGRAM ELEMENT NUMBER 

6. AUTHOR(S) 

Wilson, John W.; Slaba, Tony C.; Badavi, Francis F.; Reddell, Brandon; 
Bahadori, Amir A. 

5d. PROJECT NUMBER 

5e. TASK NUMBER 

5f. WORK UNIT NUMBER 

651549.02.07.10 

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 

NASA Langley Research Center 
Hampton, VA 23681-2199 

8. PERFORMING ORGANIZATION 
REPORT NUMBER 

L-20415 

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) 

National Aeronautics and Space Administration 
Washington, DC 20546-0001 

10. SPONSOR/MONITOR'S ACRONYM(S) 

NASA 

11. SPONSOR/MONITOR'S REPORT 
NUMBER(S) 

NASA/TP-2014-218271 


12. DISTRIBUTION/AVAILABILITY STATEMENT 

Unclassified - Unlimited 
Subject Category 93 

Availability: NASA CASI (443) 757-5802 

13. SUPPLEMENTARY NOTES 


14. ABSTRACT 

The computationally efficient HZETRN code has been used in recent trade studies for lunar and Martian exploration and is 
currently being used in the engineering development of the next generation of space vehicles, habitats, and extra vehicular 
activity equipment. A new version (3DHZETRN) capable of transporting High charge (Z) and Energy (HZE) and light ions 
(including neutrons) under space-like boundary conditions with enhanced neutron and light ion propagation is under 
development. In the present report, new algorithms for light ion and neutron propagation with well-defined convergence 
criteria in 3D objects is developed and tested against Monte Carlo simulations to verify the solution methodology. The code 
will be available through the software system, OLTARIS, for shield design and validation and provides a basis for personal 
computer software capable of space shield analysis and optimization. 


15. SUBJECT TERMS 


HZETRN; Particle transport; Radiation 


16. SECURITY CLASSIFICATION OF: 

17. LIMITATION OF 
ABSTRACT 

18. NUMBER 
OF 

19a. NAME OF RESPONSIBLE PERSON 

a. REPORT 

b. ABSTRACT 

c. THIS PAGE 

PAGES 

STI Help Desk (email: help@sti.nasa.gov) 






19b. TELEPHONE NUMBER (Include area code) 

U 

U 

u 

UU 

40 

(443) 757-5802 


Standard Form 298 (Rev. 8-98) 

Prescribed by ANSI Std. Z39.18 



































